function LaplaceJacobitest
% File name: LaplaceJacobitest.m
% Sets up and solves Laplace's equation for the test problem of Chapter 3
% using the Jacobi iterative approach.
% The problem to solve is b = A u where u contains the unknown values at the interior
% grid points.
% *** this program assumes that the grid spacing is the same in x and y directions ***
%
% Relative to the standard traverse of the unknowns in Figure 3.2:
% u = [ u(2,2)
% u(3,2)
% u(4,2)
% u(2,3)
% u(3,3)
% u(4,3) ]
%
% subfunction: fdistance
%
clear; clc   % clear variables, console and figures
%
tol=0.5*10^(-3); % tolerance for the iteration
maxiterations=50; % maximum number of iterations
NX=5; % number of grid points in x direction
NY=4; % number of grid points in y direction
N=(NX-2)*(NY-2); % number of unknowns.
%
% Create a matrix for the initial values of u where the unknown values
% are set to zero
u0=zeros(NX,NY); % current value
%
% Fill in the fixed boundary values for the problem:
u0(1,1)=6.1;
u0(2,1)=6.8;
u0(3,1)=7.7;
u0(4,1)=8.7;
u0(5,1)=9.8;
u0(1,2)=7.2;
u0(5,2)=9.4;
u0(1,3)=8.4;
u0(5,3)=9.2;
u0(1,4)=8.9;
u0(2,4)=8.9;
u0(3,4)=8.9;
u0(4,4)=8.9;
u0(5,4)=8.9;
% Create a matrix for the next iteration of u:
u1=u0; % next values
%
disp('Solving the system of linear equations for u by Jacobi iteration')
for n=1:maxiterations % iteration loop
for j=2:NY-1;
for i=2:NX-1;
u1(i,j)=(u0(i,j-1)+u0(i-1,j)+u0(i+1,j)+u0(i,j+1))/4;
end
end
d=fdistance(u1,u0); % distance between solutions at consecutive iterations
if (d < tol)
break
end
u0=u1; % copy updated solution into current solution matrix
end % of iteration loop
disp('number of iterations=')
n
disp('reshape u0 so that it looks like the data in Figure 3.2')
u0=u0'; % transpose u0 then write down its rows in reverse order
u0=u0(NY:-1:1,:);
disp('u0=')
u0
% end of LaplaceJacobitest.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d=fdistance(A,B)
% Author: Clive Mingham
% Date: 24/9/2010
% Description: Function to compute the distance between 2 N-vectors using
% the infinity norm.
[ma,na]=size(A);
[mb,nb]=size(B);
if ((ma ~= mb)|(na ~= nb) )
disp('matrices are not of equal size, stopping program')
exit
end
diff=A-B;
d=max(max(diff)); % diff is a matrix so need 2 max functions
% end of fdistance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%